


/*
This program generates alternative IV results:
Other firms' abatement decisions;


*/




*------------------------------------------------------------
*Alternative IV1: others' abatement
*----------------------------------------------------------

use  "D:\Nanjing\2019\pollution2\submission\files for submission\data\tables_comp",clear

rename exp_dummy2 Export
rename so2_new1 Abatement_so2
rename dust_new1 Abatement_dust
rename tfp_acf_new Productivity
rename lntariff Tariff
rename lnkl_lag lnKL
rename treatment_tfp tariff_prod
*rename treatment_so2 ventilation_prod
ivreghdfe lnso2  (Export Abatement_so2= tariff_prod IV2_so2 IV1_so2)   Productivity Tariff,   a(cic prov year ownership) cluster (prov_c) 
est sto so2_IV1


ivreghdfe lnso2  (Export Abatement_so2= tariff_prod IV2_so2 IV1_so2)   Productivity Tariff Tariff lnKL lnso2_ini,   a(cic prov year ownership) cluster (prov_c) 
est sto so2_IV2 


ivreghdfe lndust  (Export Abatement_so2= tariff_prod IV2_dust IV1_dust)   Productivity Tariff,   a(cic prov year ownership) cluster (prov_c) 
est sto so2_IV1

ivreghdfe lndust  (Export Abatement_dust= tariff_prod IV2_dust IV1_dust)   Productivity Tariff Tariff lnKL lndust_ini,   a(cic prov year ownership) cluster (prov_c) 
est sto dust_IV2 

 
/*
* compute R2
predict y_hat, xb
egen y_bar=mean(lnso2)

gen diff1=(y_hat-y_bar)^2
egen SSR=sum(diff1)

gen diff2=(lndust-y_bar)^2
egen SST=sum(diff2)

gen R2=SSR/SST 
sum R2
*/

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
use  "D:\Nanjing\2019\pollution2\submission\files for submission\data\tables_comp2",clear
rename exp_dummy2 Export
rename so2_new1 Abatement_so2
rename dust_new1 Abatement_dust
rename tfp_acf_new Productivity
rename lntariff Tariff
rename lnkl_lag lnKL
rename treatment_tfp tariff_prod
rename treatment_so2 ventilation_prod




 capture drop ub
 egen ub=pctile(lndust),p(97)
 drop if lnKL==-2.5

 ***********************************************************
*  Columns 3 and 7 of Table 4 (first stage in Table 3)
*********************************************************


*egen prov_d=group(party_id year cic indc)
 
 ivreghdfe lnso2  (Export Abatement_so2= tariff_prod IV2_so2 IV1_so2)  ///
 Productivity Tariff lnKL  , a(indc year ) cluster(prov_c) 

 
 ivreghdfe lndust  (Export Abatement_dust= tariff_prod IV2_dust IV1_dust)  Productivity Tariff lnKL if lndust<ub,  a(indc  year) cluster(prov_c) 


 
 
drop Productivity Tariff lnKL Abatement_so2 Abatement_dust
 
rename so2_new1_d Abatement_so2
rename dust_new1_d Abatement_dust
rename tfp_acf_new_d Productivity
rename lntariff_d Tariff
rename lnkl_d lnKL


ivreghdfe lnso2_d  (Export Abatement_so2=  IV1_so2 IV2_so2 tariff_prod )   Productivity Tariff lnKL , a(cic prov year ownership ) cluster(prov_c) 

ivreghdfe lndust_d  (Export Abatement_dust= IV1_dust IV2_dust   tariff_prod)  Productivity Tariff lnKL , a(cic prov year ownership) cluster(prov_c) 
***************
 /*
* compute R2
predict y_hat, xb
egen y_bar=mean(lnso2)

gen diff1=(y_hat-y_bar)^2
egen SSR=sum(diff1)

gen diff2=(lnso2_d-y_bar)^2
egen SST=sum(diff2)

gen R2=SSR/SST 
sum R2
*/



















 
 
 
 
 *--------------------------------------------------------------------------------
 *--------------------------------------------
 *Bootstrap for alternative IV
 *---------------------------------------------
 
 
 

 clear mata
   mata:

void GMM_DLW_CD(todo,betas,crit,g,H)
  {
	PHI=st_data(.,("phi"))
	PHI_LAG=st_data(.,("phi_lag"))
	Z=st_data(.,("const","m_lag","l_lag","y1","y2","y3","k","k_lag","exp_dummy99_lag","exp_dummy00_lag","exp_dummy01_lag","exp_dummy02_lag","exp_dummy03_lag","exp_dummy04_lag","exp_dummy05_lag"))
	X=st_data(.,("const","m","l","x1","x2","x3","k","exp_dummy99","exp_dummy00","exp_dummy01","exp_dummy02","exp_dummy03","exp_dummy04","exp_dummy05"))
	X_lag=st_data(.,("const","m_lag","l_lag","x1_lag","x2_lag","x3_lag","k_lag","exp_dummy99_lag","exp_dummy00_lag","exp_dummy01_lag","exp_dummy02_lag","exp_dummy03_lag","exp_dummy04_lag","exp_dummy05_lag"))
	Y=st_data(.,("q"))
	C=st_data(.,("const"))
	

	EXPORT=st_data(.,("exp_dummy99_lag","exp_dummy00_lag","exp_dummy01_lag","exp_dummy02_lag","exp_dummy03_lag","exp_dummy04_lag","exp_dummy05_lag"))

	OMEGA=PHI-X*betas'
	OMEGA_lag=PHI_LAG-X_lag*betas'
	OMEGA_lag_pol=(C,EXPORT,OMEGA_lag)
	g_b = invsym(OMEGA_lag_pol'OMEGA_lag_pol)*OMEGA_lag_pol'OMEGA 
	XI=OMEGA-OMEGA_lag_pol*g_b /*XI=epsilon: recover epsilon_it in equation (23): shock in productivity;*/
	crit=(Z'XI)'(Z'XI)
  }
///
///
void DLW_TRANSLOG()
	{
	initialvalue=st_data(1,("initialConst","initialm","initiall","initialx1","initialx2","initialx3","initialk","initialexp_dummy99","initialexp_dummy00","initialexp_dummy01","initialexp_dummy02","initialexp_dummy03","initialexp_dummy04","initialexp_dummy05"))
	S=optimize_init()
	optimize_init_evaluator(S, &GMM_DLW_CD())
	optimize_init_evaluatortype(S,"d0")
	optimize_init_technique(S, "nm")
	optimize_init_nmsimplexdeltas(S, 0.1)
	optimize_init_which(S,"min")
	optimize_init_params(S,initialvalue)
	p=optimize(S)
	p
	st_matrix("beta_dlwtranslog",p)
	}
end

capture program drop dlw_translog
program dlw_translog, rclass
preserve
**********************************************
sort firm year
mata DLW_TRANSLOG()
end 
 use "D:\Nanjing\2019\pollution2\submission\files for submission\data\tables_comp",clear

 gen cic_adj=cic

gen hangye=int(cic_adj/100)
drop if hangye<13
drop if hangye>42

*keep if hangye==27|hangye==14
capture drop cic2
capture drop cic4 
gen cic2=floor(cic_adj/100)
egen gcic4=group(cic_adj)
*egen SDcic=sd(cic2), by(firm)
*egen cic_y=group(cic year)
gen Productivity=tfp_acf_new
*----------------------------------------------------------------------
* bootstrap starts from here
*----------------------------------------------------------------------

capture program drop prodboot
program define prodboot, rclass
preserve






xi:heckman ratio1 dyear*, noconstant twostep sel(exp_dummy=  tfp_acf_new clean_dummy_so2_lag clean_dummy_dust_lag clean_dummy_cod_lag dcic2* lnv_lag lnkl_lag dprov* owner* dyear*)  

gen coef1_bs=_b[dyear1]
gen coef2_bs=_b[dyear2]
gen coef3_bs=_b[dyear3]
gen coef4_bs=_b[dyear4]
gen coef5_bs=_b[dyear5]
gen coef6_bs=_b[dyear6]
gen coef7_bs=_b[dyear7]
capture drop exp_dummy2
gen exp_dummy2=.
replace exp_dummy2=exp_dummy*coef1_bs*dyear1 if year==1999
replace exp_dummy2=exp_dummy*coef2_bs*dyear2 if year==2000
replace exp_dummy2=exp_dummy*coef3_bs*dyear3 if year==2001
replace exp_dummy2=exp_dummy*coef4_bs*dyear4 if year==2002
replace exp_dummy2=exp_dummy*coef5_bs*dyear5 if year==2003
replace exp_dummy2=exp_dummy*coef6_bs*dyear6 if year==2004
replace exp_dummy2=exp_dummy*coef7_bs*dyear7 if year==2005
drop  coef1_bs coef2_bs coef3_bs coef4_bs coef5_bs coef6_bs coef7_bs






*************************************************************************************************
* For so2+ IV: column 1 
ivreg2 lnso2  (exp_dummy2 so2_new1= IV1_so2 IV2_so2 treatment_tfp)   ///
tfp_acf_new lntariff  i.cic i.prov i.year, partial(i.cic i.prov i.year) cluster (prov_c)

return scalar exp_t2c1=_b[exp_dummy2]
return scalar abt_t2c1=_b[so2_new1]
return scalar tfp_t2c1=_b[tfp_acf_new]
return scalar tar_t2c1=_b[lntariff]


********************************************************************************************************

* column 2
ivreg2 lnso2  (exp_dummy2 so2_new1= IV1_so2 IV2_so2 treatment_tfp)   ///
tfp_acf_new lntariff lnkl_lag  lnso2_ini i.cic i.prov i.year, partial(i.cic i.prov i.year) cluster (prov_c)

return scalar exp_t2c2=_b[exp_dummy2]
return scalar abt_t2c2=_b[so2_new1]
return scalar tfp_t2c2=_b[tfp_acf_new]
return scalar tar_t2c2=_b[lntariff]
return scalar klr_t2c2=_b[lnkl_lag]
return scalar iei_t2c2=_b[lnso2_ini]


************************************************************************************************************************************
* For dust +IV: column 5
ivreg2 lndust  (exp_dummy2 dust_new1= IV1_dust IV2_dust treatment_tfp) ///
  tfp_acf_new lntariff  i.cic i.prov i.year, partial(i.cic i.prov i.year ) cluster (prov_c)

return scalar exp_t2c3=_b[exp_dummy2]
return scalar abt_t2c3=_b[dust_new1]
return scalar tfp_t2c3=_b[tfp_acf_new]
return scalar tar_t2c3=_b[lntariff]

***************************************************************************************************************************************



************************************************************************************************************************************
* For dust +IV: column 6
ivreg2 lndust  (exp_dummy2 dust_new1= IV1_dust IV2_dust treatment_tfp) ///
  tfp_acf_new lntariff lnkl_lag  lndust_ini i.cic i.prov i.year, partial(i.cic i.prov i.year ) cluster (prov_c)

return scalar exp_t2c4=_b[exp_dummy2]
return scalar abt_t2c4=_b[dust_new1]
return scalar tfp_t2c4=_b[tfp_acf_new]
return scalar tar_t2c4=_b[lntariff]
return scalar klr_t2c4=_b[lnkl_lag]
return scalar iei_t2c4=_b[lndust_ini]
***************************************************************************************************************************************


end

capture drop indc
egen indc=group(name_raw)
drop if indc==.

*-------------------------------------------------------
* For columns 1, 2, 5, 6
*--------------------------------------------------------

bootstrap  exp_t2c1=r(exp_t2c1) abt_t2c1=r(abt_t2c1) tfp_t2c1=r(tfp_t2c1) tar_t2c1=r(tar_t2c1)  ///
 exp_t2c2=r(exp_t2c2) abt_t2c2=r(abt_t2c2) tfp_t2c2=r(tfp_t2c2) tar_t2c2=r(tar_t2c2)  klr_t2c2=r(klr_t2c2) iei_t2c2=r(iei_t2c2)    ///
 exp_t2c3=r(exp_t2c3) abt_t2c3=r(abt_t2c3) tfp_t2c3=r(tfp_t2c3) tar_t2c3=r(tar_t2c3)   ///
 exp_t2c4=r(exp_t2c4) abt_t2c4=r(abt_t2c4) tfp_t2c4=r(tfp_t2c4) tar_t2c4=r(tar_t2c4)  klr_t2c4=r(klr_t2c4) iei_t2c4=r(iei_t2c2)   ,  ///
 reps(50)  seed(1) nodrop cluster(indc) idcluster(newid): prodboot

 
 
 
 ********************************************************************************
 
 
 * for column 3,4, 7, 8
 

 * for column 3,4, 7, 8
 
  clear mata
   mata:

void GMM_DLW_CD(todo,betas,crit,g,H)
  {
	PHI=st_data(.,("phi"))
	PHI_LAG=st_data(.,("phi_lag"))
	Z=st_data(.,("const","m_lag","l_lag","y1","y2","y3","k","k_lag","exp_dummy99_lag","exp_dummy00_lag","exp_dummy01_lag","exp_dummy02_lag","exp_dummy03_lag","exp_dummy04_lag","exp_dummy05_lag"))
	X=st_data(.,("const","m","l","x1","x2","x3","k","exp_dummy99","exp_dummy00","exp_dummy01","exp_dummy02","exp_dummy03","exp_dummy04","exp_dummy05"))
	X_lag=st_data(.,("const","m_lag","l_lag","x1_lag","x2_lag","x3_lag","k_lag","exp_dummy99_lag","exp_dummy00_lag","exp_dummy01_lag","exp_dummy02_lag","exp_dummy03_lag","exp_dummy04_lag","exp_dummy05_lag"))
	Y=st_data(.,("q"))
	C=st_data(.,("const"))
	

	EXPORT=st_data(.,("exp_dummy99_lag","exp_dummy00_lag","exp_dummy01_lag","exp_dummy02_lag","exp_dummy03_lag","exp_dummy04_lag","exp_dummy05_lag"))

	OMEGA=PHI-X*betas'
	OMEGA_lag=PHI_LAG-X_lag*betas'
	OMEGA_lag_pol=(C,EXPORT,OMEGA_lag)
	g_b = invsym(OMEGA_lag_pol'OMEGA_lag_pol)*OMEGA_lag_pol'OMEGA 
	XI=OMEGA-OMEGA_lag_pol*g_b /*XI=epsilon: recover epsilon_it in equation (23): shock in productivity;*/
	crit=(Z'XI)'(Z'XI)
  }
///
///
void DLW_TRANSLOG()
	{
	initialvalue=st_data(1,("initialConst","initialm","initiall","initialx1","initialx2","initialx3","initialk","initialexp_dummy99","initialexp_dummy00","initialexp_dummy01","initialexp_dummy02","initialexp_dummy03","initialexp_dummy04","initialexp_dummy05"))
	S=optimize_init()
	optimize_init_evaluator(S, &GMM_DLW_CD())
	optimize_init_evaluatortype(S,"d0")
	optimize_init_technique(S, "nm")
	optimize_init_nmsimplexdeltas(S, 0.1)
	optimize_init_which(S,"min")
	optimize_init_params(S,initialvalue)
	p=optimize(S)
	p
	st_matrix("beta_dlwtranslog",p)
	}
end

capture program drop dlw_translog
program dlw_translog, rclass
preserve
**********************************************
sort firm year
mata DLW_TRANSLOG()
end 
use "D:\Nanjing\2019\pollution2\submission\files for submission\data\tables_comp2",clear
gen Export=exp_dummy2 
* gen frdm=party_id
 gen cic_adj=cic
*rename frdm firm
gen hangye=int(cic_adj/100)
drop if hangye<13
*drop if hangye<40
drop if hangye>42

*keep if hangye==27|hangye==14
capture drop cic2
capture drop cic4 
gen cic2=floor(cic_adj/100)
egen gcic4=group(cic_adj)
*egen SDcic=sd(cic2), by(firm)
capture drop cic_y
egen cic_y=group(cic year)
capture drop k



gen indc_firm=indc
gen indc2 =indc_firm

*rename party_id frdm

*
*----------------------------------------------------------------------
* bootstrap starts from here
*----------------------------------------------------------------------
*keep if _merge==2|_merge==3
capture program drop prodboot
program define prodboot, rclass
preserve






* Re-estimate firm-level production function for each sample from bootstrap



xi:heckman ratio1 dyear*, noconstant twostep sel(exp_dummy=  tfp_acf_new clean_dummy_so2_lag clean_dummy_dust_lag clean_dummy_cod_lag dcic2* lnv_lag lnkl_lag dprov* owner* dyear*)  

gen coef1_bs=_b[dyear1]
gen coef2_bs=_b[dyear2]
gen coef3_bs=_b[dyear3]
gen coef4_bs=_b[dyear4]
gen coef5_bs=_b[dyear5]
gen coef6_bs=_b[dyear6]
gen coef7_bs=_b[dyear7]
capture drop exp_dummy2
gen exp_dummy2=.
replace exp_dummy2=exp_dummy*coef1_bs*dyear1 if year==1999
replace exp_dummy2=exp_dummy*coef2_bs*dyear2 if year==2000
replace exp_dummy2=exp_dummy*coef3_bs*dyear3 if year==2001
replace exp_dummy2=exp_dummy*coef4_bs*dyear4 if year==2002
replace exp_dummy2=exp_dummy*coef5_bs*dyear5 if year==2003
replace exp_dummy2=exp_dummy*coef6_bs*dyear6 if year==2004
replace exp_dummy2=exp_dummy*coef7_bs*dyear7 if year==2005
drop  coef1_bs coef2_bs coef3_bs coef4_bs coef5_bs coef6_bs coef7_bs







gen Abatement_so2=so2_new1 
gen Abatement_dust= dust_new1 
gen  Productivity= tfp_acf_new
gen Tariff= lntariff 
gen lnKL= lnkl_lag 
gen tariff_prod=treatment_tfp 

*************************************************************************************************
* For so2+ OL: column 3 of Table 2
ivreghdfe lnso2  (Export Abatement_so2= tariff_prod IV1_so2 IV2_so2) Productivity Tariff lnKL , a(indc_firm year)

return scalar exp_t1c2=_b[Export]
return scalar abt_t1c2=_b[Abatement_so2]
return scalar tfp_t1c2=_b[Productivity]
return scalar tar_t1c2=_b[Tariff]
return scalar klr_t1c2=_b[lnKL]

capture drop ub
egen ub=pctile(lndust), p(97)
 ivreghdfe lndust  (Export Abatement_dust= tariff_prod IV1_dust IV2_dust) Productivity Tariff lnKL if lndust<ub, a(indc_firm year)

return scalar exp_t2c2=_b[Export]
return scalar abt_t2c2=_b[Abatement_dust]
return scalar tfp_t2c2=_b[Productivity]
return scalar tar_t2c2=_b[Tariff]
return scalar klr_t2c2=_b[lnKL]


 
 
drop Productivity Tariff lnKL Abatement_so2 Abatement_dust
 
rename so2_new1_d Abatement_so2
rename dust_new1_d Abatement_dust
rename tfp_acf_new_d Productivity
rename lntariff_d Tariff
rename lnkl_d lnKL

ivreghdfe lnso2_d  (Export Abatement_so2=  IV1_so2 IV2_so2 tariff_prod )   Productivity Tariff lnKL , a(cic prov year ownership ) cluster(prov_c) 

return scalar exp_t3c2=_b[Export]
return scalar abt_t3c2=_b[Abatement_so2]
return scalar tfp_t3c2=_b[Productivity]
return scalar tar_t3c2=_b[Tariff]
return scalar klr_t3c2=_b[lnKL]

ivreghdfe lndust_d  (Export Abatement_dust= IV1_dust IV2_dust   IV1_so2)  Productivity Tariff lnKL , a(cic prov year ownership) cluster(prov_c) 

return scalar exp_t4c2=_b[Export]
return scalar abt_t4c2=_b[Abatement_dust]
return scalar tfp_t4c2=_b[Productivity]
return scalar tar_t4c2=_b[Tariff]
return scalar klr_t4c2=_b[lnKL]
***************
end

*-------------------------------------------------------
* For columns 3, 4, 7, 8
*--------------------------------------------------------

bootstrap   exp_t1c2=r(exp_t1c2) abt_t1c2=r(abt_t1c2) tfp_t1c2=r(tfp_t1c2) tar_t1c2=r(tar_t1c2)  klr_t1c2=r(klr_t1c2) ///
exp_t2c2=r(exp_t2c2) abt_t2c2=r(abt_t2c2) tfp_t2c2=r(tfp_t2c2) tar_t2c2=r(tar_t2c2)  klr_t2c2=r(klr_t2c2) ///
exp_t3c2=r(exp_t3c2) abt_t3c2=r(abt_t3c2) tfp_t3c2=r(tfp_t3c2) tar_t3c2=r(tar_t3c2)  klr_t3c2=r(klr_t3c2) ///
exp_t4c2=r(exp_t4c2) abt_t4c2=r(abt_t4c2) tfp_t4c2=r(tfp_t4c2) tar_t4c2=r(tar_t4c2)  klr_t4c2=r(klr_t4c2), ///
 rep(50)  seed(5) nodrop : prodboot

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
